home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / umich / falcon / programm.ing / nt_dsp1.lzh / NT_DSP1.MSA / FFT / DCT1.ASM next >
Assembly Source File  |  1990-01-17  |  6KB  |  149 lines
  1. ;
  2. ; This program, originally available on the Motorola DSP bulletin board,
  3. ; is provided under a DISCLAIMER OF WARRANTY available from Motorola DSP
  4. ; Operation, 6501 William Cannon Drive, West, Austin, Texas  78735-8598.
  5. ;
  6.         ident   1,1
  7.         page    132,56,1,1
  8.  
  9.         opt     mex,nomd,cex
  10. ;
  11. ;       Discrete Cosine Transform (DCT)
  12. ;
  13. ;    V1.2
  14. ;    03 Oct 88
  15. ;
  16. ;       The DCT is implemented by using the method described in
  17. ;       "On the Computation of the Discrete Cosine Transform", IEEE
  18. ;       Transactions on Communications. Vol COM-26, No. 6, June 1978
  19. ;
  20. ;       The operation of this transform is as follows:
  21. ;           1.  Resort the original sequence.
  22. ;           2.  Perform an FFT on the modified sequence.
  23. ;           3.  Multiply the output of the FFT by a complex exponential
  24. ;               and take the real part.
  25. ;           4.  Perform a bit reversed sort to unscramble the result.
  26. ;
  27. ;       After the transform, the data needs to be scaled (according to
  28. ;       the previously cited reference):
  29. ;           g(0) = g(0)*1.414/N
  30. ;           g(k) = g(k)*2.0/N       k=1,2,...,N-1
  31. ;       This scaling is NOT performed by this program.  This type of scaling
  32. ;       is application dependent and some implementations may not require
  33. ;       the scaling.
  34. ;
  35.  
  36.  
  37.         include 'sincos'
  38.         include 'fftr2aa'
  39.  
  40. start   equ     $100            ;beginning of program in P memory
  41. points  equ     16              ;number of points in dct
  42. data    equ     0               ;address of data in Y memory
  43. coef    equ     16              ;address of coefficients
  44.  
  45.         sincos  points,coef     ;generate twiddle factor for FFT.
  46.  
  47. ;       Generate exponential table after twiddle factor table.
  48. ;
  49. ;       This generates the scale factor of EXP(-j*pi*k/(2*N)).
  50. ;       The real part is in Y memory and the imaginary part is
  51. ;       in X memory.  Note that the negative of the part in Y memory
  52. ;       is actually stored so that a -1 (exactly) can be represented
  53. ;       with the fractional arithmetic.  This negative sign is
  54. ;       compensated when the actual rotation is performed.
  55. ;
  56. rot     equ     pi/@cvf(2*points)
  57.  
  58.         org     x:
  59. exptbl
  60. count   set     0
  61.         dup     points
  62.         dc      -@sin(@cvf(count)*rot)
  63. count   set     count+1
  64.         endm
  65.  
  66.         org     y:
  67. count   set     0
  68.         dup     points
  69.         dc      -@cos(@cvf(count)*rot)
  70. count   set     count+1
  71.         endm
  72.  
  73.  
  74.  
  75.         org     p:start
  76. ;
  77. ;       Do inital sort. The input data is in Y, the sorted data goes to X.
  78. ;       This forms a new sequence in X memory from the sequence in Y memory
  79. ;       according to the definition:
  80. ;               x(k)=y(2k)
  81. ;               x(N-1-K)=x(2k+1)      k=0,1,...,N/2-1
  82. ;
  83.         move    #$ffff,m0               ;set linear
  84.         move    m0,m4                   ;ditto
  85.         move    #data,r0                ;point to data
  86.         move    #2,n0                   ;jump by even numbers
  87.         move    r0,r4                   ;copy data pointer
  88.         move    y:(r0)+n0,a             ;get first point
  89.         rep     #points/2               ;sort even points
  90.         move    y:(r0)+n0,a  a,x:(r4)+
  91.         move    #data+1,r0              ;point to odds
  92.         move    #data+points-1,r4       ;odds output
  93.         move    y:(r0)+n0,a             ;get first odd
  94.         rep     #points/2               ;sort odd points
  95.         move    y:(r0)+n0,a  a,x:(r4)-
  96.         clr     a       #data,r0        ;get a zero, point r0 to data
  97.         rep     #points                 ;clear out original data
  98.         move    a,y:(r0)+
  99.  
  100. ;
  101. ;       Do an FFT on the new sequence
  102. ;
  103.         fftr2aa  points,data,coef        ;do FFT
  104.  
  105. ;
  106. ;       Multiply the output by a complex exponential and take the
  107. ;       real part.  dct=re[ (a+jb)*(c+jd)] = a*c-b*d where a+jb is the
  108. ;       output of the FFT and c+jd is the rotation factor. Note that
  109. ;       the data is still in bit reversed order but the rotation
  110. ;       table is in sequential order.  Thus, the data is accessed with
  111. ;       bit reversed addressing but the rotation table is addressed
  112. ;       with linear addressing.
  113. ;
  114.         move    #exptbl,r0              ;point to exponential table
  115.         move    #$ffff,m0               ;linear addressing
  116.  
  117.         move    #data,r4                ;point to output of fft
  118.         move    r4,r5                   ;copy data pointer
  119.         move    #0,m4                   ;bit reversed addressing
  120.         move    m4,m5
  121.         move    #points/2,n4
  122.         move    n4,n5
  123.  
  124.         move    x:(r4),x0  y:(r0),y0                    ;get a,-c
  125.         do      #points,_gf                             ;do all points
  126.         mpy     -x0,y0,a  x:(r0)+,x1   y:(r4)+n4,y1     ;-a*c, get d,b
  127.         macr    -x1,y1,a  x:(r4),x0    y:(r0),y0        ;-b*d, get a,c
  128.         move              a,x:(r5)+n5                   ;save dct
  129. _gf
  130.  
  131. ;
  132. ;       Bit reverse the output. The rotated and bit reversed data from X
  133. ;       memory is transfered to Y memory and unscrambled.  The resulting
  134. ;       DCT is in Y memory starting at DATA with length POINTS.
  135. ;
  136.         move    #data,r0                ;point to data
  137.         move    #points/2,n0            ;number of points/2
  138.         move    #0,m0                   ;set reverse carry
  139.  
  140.         move    r0,r4                   ;copy pointer
  141.         move    #$ffff,m4               ;set linear
  142.  
  143.         move    x:(r0)+n0,a             ;do preload
  144.         rep     #points                 ;unscramble
  145.         move    x:(r0)+n0,a   a,y:(r4)+
  146.  
  147.         end
  148. ^Z
  149.